####Model results: adding interaction autonomy x included (civil war) (table X7, models A22-A23)####
main_m2g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c(m2g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
main_m2g_t1_sc_cse <- data.frame(cluster.se(main_m2g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
main_m2g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c(m2g_t1, group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
main_m2g_t1_sf_cse <- data.frame(cluster.se(main_m2g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
stargazer(main_m2g_t1_sc, main_m2g_t1_sf, se=c(main_m2g_t1_sc_cse, main_m2g_t1_sf_cse), dep.var.labels.include = F, omit=c("cowcode|region|factor|years"), type="text", order=vars.order_g, title="Table X7. Additional results: interacting territorial autonomy and included.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil violence involving a group in a given unit.","Region- and year-fixed effects included but not reported."), model.numbers =F, column.labels = c("Maj. (model A22)", "Min. (model A23)"), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included","Included","Relative size (state)","Relative size (unit)", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil war spatial lag"), out="../tables/additional_results_tables_appendices/tablex7.txt")
stargazer(main_m2g_t1_sc, main_m2g_t1_sf, se=c(main_m2g_t1_sc_cse, main_m2g_t1_sf_cse), dep.var.labels.include = F, omit=c("cowcode|region|factor|years"), type="html", order=vars.order_g, title="Table X7. Additional results: interacting territorial autonomy and included.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil violence involving a group in a given unit.","Region- and year-fixed effects included but not reported."), model.numbers =F, column.labels = c("Maj. (model A22)", "Min. (model A23)"), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included","Included","Relative size (state)","Relative size (unit)", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil war spatial lag"), out="../tables/additional_results_tables_appendices/tablex7.html")

####Model results: distinguishing whether SO maj or min is excluded / distinguishing powerless/discriminated / both (table X8, models A24-A27)####
main_m3d_t1_sb <- glm(as.formula(paste("cv_event", paste(c(m3d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
main_m3d_t1_sb_cse <- data.frame(cluster.se(main_m3d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r6_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t","included_powerless","included_discriminated","powerless_powerless","powerless_discriminated","discriminated_discriminated", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r6_m1d_t1_sb_cse <- data.frame(cluster.se(r6_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r6_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t*included_powerless","sa_territory_t*included_discriminated","powerless_powerless","powerless_discriminated","discriminated_discriminated", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r6_m2d_t1_sb_cse <- data.frame(cluster.se(r6_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r6_m3d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t*includedd_powerlessnd","sa_territory_t*includednd_powerlessd","sa_territory_t*includedd_discriminatednd","sa_territory_t*includednd_discriminatedd","powerless_powerless","powerless_discriminated","discriminated_discriminated", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r6_m3d_t1_sb_cse <- data.frame(cluster.se(r6_m3d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(main_m3d_t1_sb, r6_m1d_t1_sb, r6_m2d_t1_sb, r6_m3d_t1_sb, se=c(main_m3d_t1_sb_cse, r6_m1d_t1_sb_cse, r6_m2d_t1_sb_cse, r6_m3d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Communal violence","Communal violence","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj./min. dyad (model A24)","Maj./min. dyad (model A25)","Maj./min. dyad (model A26)","Maj./min. dyad (model A27)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X8. Additional results: identifying whether second-order minority or majority is included/excluded; distinguishing between powerless and discriminated exclusion.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of communal violence involving a dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included (maj.)/excluded (min.)","Territorial autonomy x included (min.)/excluded (maj.)","Territorial autonomy x included/powerless","Territorial autonomy x included/discriminated","Territorial autonomy x included (maj.)/powerless (min.)","Territorial autonomy x included (min.)/powerless (maj.)","Territorial autonomy x included (maj.)/discriminated (min.)","Territorial autonomy x included (min.)/discriminated (maj.)","Included (maj.)/excluded (min.)","Included (min.)/excluded (maj.)","Included/powerless","Included/discriminated","Included (maj.)/powerless (min.)","Included (min.)/powerless (maj.)","Included (maj.)/discriminated (min.)","Included (min.)/discriminated (maj.)","Excluded/excluded","Powerless/powerless","Powerless/discriminated","Discriminated/discriminated","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex8.txt")
stargazer(main_m3d_t1_sb, r6_m1d_t1_sb, r6_m2d_t1_sb, r6_m3d_t1_sb, se=c(main_m3d_t1_sb_cse, r6_m1d_t1_sb_cse, r6_m2d_t1_sb_cse, r6_m3d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Communal violence","Communal violence","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj./min. dyad (model A24)","Maj./min. dyad (model A25)","Maj./min. dyad (model A26)","Maj./min. dyad (model A27)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X8. Additional results: identifying whether second-order minority or majority is included/excluded; distinguishing between powerless and discriminated exclusion.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of communal violence involving a dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included (maj.)/excluded (min.)","Territorial autonomy x included (min.)/excluded (maj.)","Territorial autonomy x included/powerless","Territorial autonomy x included/discriminated","Territorial autonomy x included (maj.)/powerless (min.)","Territorial autonomy x included (min.)/powerless (maj.)","Territorial autonomy x included (maj.)/discriminated (min.)","Territorial autonomy x included (min.)/discriminated (maj.)","Included (maj.)/excluded (min.)","Included (min.)/excluded (maj.)","Included/powerless","Included/discriminated","Included (maj.)/powerless (min.)","Included (min.)/powerless (maj.)","Included (maj.)/discriminated (min.)","Included (min.)/discriminated (maj.)","Excluded/excluded","Powerless/powerless","Powerless/discriminated","Discriminated/discriminated","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex8.html")

####Figure A7####
##a) adding interaction autonomy x included (civil war)
me_main_m2g_t1_sc <- margins_summary(main_m2g_t1_sc, variables = "sa_territory_t", at = list(included = c(0,1)), vcov = cluster.vcov(main_m2g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_main_m2g_t1_sc$term <- ifelse(me_main_m2g_t1_sc$included == 1, "maj., included", "maj., excluded")
me_main_m2g_t1_sf <- margins_summary(main_m2g_t1_sf, variables = "sa_territory_t", at = list(included = c(0,1)), vcov = cluster.vcov(main_m2g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_main_m2g_t1_sf$term <- ifelse(me_main_m2g_t1_sf$included == 1, "min., included", "min., excluded")
me_main_m2g <- rbind.fill(me_main_m2g_t1_sc, me_main_m2g_t1_sf)
me_main_m2g$term <- as.factor(me_main_m2g$term)
me_main_m2g$term = factor(me_main_m2g$term,levels(me_main_m2g$term)[c(4,3,2,1)])
me_main_m2g$model <- 1
me_main_m2g$estimate <- me_main_m2g$AME
me_main_m2g$conf.low <- me_main_m2g$lower
me_main_m2g$conf.high <- me_main_m2g$upper
me_main_m2g_plot <- dwplot(me_main_m2g, vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2), dot_args = list(shape = 16, colour = "black"), whisker_args = list(linetype = 1, colour = "black")) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence (TA x inclusion)") +
  scale_color_manual(name = "Coefficient for:", values = c("#1b9e77","#d95f02","#7570b3"))+ scale_linetype_manual(name = "Coefficient for:", values = c(3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+
  guides(shape = guide_legend(nrow=2,"model",reverse=TRUE), colour = guide_legend(nrow=2,"model",reverse=TRUE),linetype = guide_legend(nrow=2,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format(accuracy = 0.5),breaks=c(-0.03,-0.015,0))
##b) distinguishing whether SO maj or min is excluded
me_main_m3d_t1_sb <- margins_summary(main_m3d_t1_sb, variables = "sa_territory_t", at = list(includedd_excludednd = c(0,1), includednd_excludedd = c(0,1)), vcov = cluster.vcov(main_m3d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_main_m3d_t1_sb <- subset(me_main_m3d_t1_sb, !(includedd_excludednd == 1 & includednd_excludedd == 1))
me_main_m3d_t1_sb$term <- ifelse(me_main_m3d_t1_sb$includedd_excludednd == 1, "maj. incl./min. excl.", NA)
me_main_m3d_t1_sb$term <- ifelse(is.na(me_main_m3d_t1_sb$term) & me_main_m3d_t1_sb$includednd_excludedd == 1, "min. incl./maj. excl.", me_main_m3d_t1_sb$term)
me_main_m3d_t1_sb$term <- ifelse(is.na(me_main_m3d_t1_sb$term), "other", me_main_m3d_t1_sb$term)
me_main_m3d_t1_sb$model <- 5
me_main_m3d_t1_sb <- subset(me_main_m3d_t1_sb, term != "other")
me_main_m3d_t1_sb$term <- as.factor(me_main_m3d_t1_sb$term)
me_main_m3d_t1_sb$term = factor(me_main_m3d_t1_sb$term,levels(me_main_m3d_t1_sb$term)[c(2,1)])
me_main_m3d_t1_sb <- me_main_m3d_t1_sb%>% arrange(desc(term))
me_main_m3d_t1_sb$estimate <- me_main_m3d_t1_sb$AME
me_main_m3d_t1_sb$conf.low <- me_main_m3d_t1_sb$lower
me_main_m3d_t1_sb$conf.high <- me_main_m3d_t1_sb$upper
me_main_m3d_t1_sb_plot <- dwplot(me_main_m3d_t1_sb,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(shape = 16, colour = "black"), whisker_args = list(linetype = 1, colour = "black")) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) + ggtitle("b) communal violence maj./min. dyad\n(maj./min. incl./excl.)") +
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type") + scale_x_continuous(labels=scales::percent_format())
##c) distinguishing powerless/discriminated / both
me_r6_m1d_t1_sb <- margins_summary(r6_m1d_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(r6_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m1d_t1_sb$term <- "all"
me_r6_m2d_t1_sb1 <- margins_summary(r6_m2d_t1_sb, variables = "sa_territory_t", at = list(included_powerless = 1), vcov = cluster.vcov(r6_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m2d_t1_sb2 <- margins_summary(r6_m2d_t1_sb, variables = "sa_territory_t", at = list(included_discriminated = 1), vcov = cluster.vcov(r6_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m2d_t1_sb <- rbind.fill(me_r6_m2d_t1_sb1,me_r6_m2d_t1_sb2)
me_r6_m2d_t1_sb$term <- ifelse(me_r6_m2d_t1_sb$included_powerless == 1, "included/powerless", NA)
me_r6_m2d_t1_sb$term <- ifelse(is.na(me_r6_m2d_t1_sb$term) & me_r6_m2d_t1_sb$included_discriminated == 1, "included/discriminated",me_r6_m2d_t1_sb$term)
me_r6_m3d_t1_sb1 <- margins_summary(r6_m3d_t1_sb, variables = "sa_territory_t", at = list(includedd_powerlessnd = 1), vcov = cluster.vcov(r6_m3d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m3d_t1_sb2 <- margins_summary(r6_m3d_t1_sb, variables = "sa_territory_t", at = list(includednd_powerlessd = 1), vcov = cluster.vcov(r6_m3d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m3d_t1_sb3 <- margins_summary(r6_m3d_t1_sb, variables = "sa_territory_t", at = list(includedd_discriminatednd = 1), vcov = cluster.vcov(r6_m3d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m3d_t1_sb4 <- margins_summary(r6_m3d_t1_sb, variables = "sa_territory_t", at = list(includednd_discriminatedd = 1), vcov = cluster.vcov(r6_m3d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r6_m3d_t1_sb <- rbind.fill(me_r6_m3d_t1_sb1,me_r6_m3d_t1_sb2,me_r6_m3d_t1_sb3,me_r6_m3d_t1_sb4)
me_r6_m3d_t1_sb$term <- ifelse(me_r6_m3d_t1_sb$includedd_powerlessnd == 1, "maj. incl./min. powerless", NA)
me_r6_m3d_t1_sb$term <- ifelse(is.na(me_r6_m3d_t1_sb$term) & me_r6_m3d_t1_sb$includednd_powerlessd == 1, "min. incl./maj. powerless",me_r6_m3d_t1_sb$term)
me_r6_m3d_t1_sb$term <- ifelse(is.na(me_r6_m3d_t1_sb$term) & me_r6_m3d_t1_sb$includedd_discriminatednd == 1, "maj. incl./min. discriminated",me_r6_m3d_t1_sb$term)
me_r6_m3d_t1_sb$term <- ifelse(is.na(me_r6_m3d_t1_sb$term) & me_r6_m3d_t1_sb$includednd_discriminatedd == 1, "min. incl./maj. discriminated",me_r6_m3d_t1_sb$term)
me_r6_m13d <- rbind.fill(me_r6_m1d_t1_sb, me_r6_m2d_t1_sb, me_r6_m3d_t1_sb)
me_r6_m13d <- subset(me_r6_m13d, term != "other")
me_r6_m13d$term <- as.factor(me_r6_m13d$term)
me_r6_m13d$term = factor(me_r6_m13d$term,levels(me_r6_m13d$term)[c(1,3,2,5,7,4,6)])
me_r6_m13d <- me_r6_m13d %>% arrange(factor(term))
me_r6_m13d$estimate <- me_r6_m13d$AME
me_r6_m13d$conf.low <- me_r6_m13d$lower
me_r6_m13d$conf.high <- me_r6_m13d$upper
me_r6_m13d_plot <- dwplot(me_r6_m13d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(shape = 16, colour = "black"), whisker_args = list(linetype = 1, colour = "black")) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("c) communal violence maj./min. dyad\n(powerless/discriminated)") +
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type") +
  scale_x_continuous(labels=scales::percent_format(),breaks=c(0,0.0075,0.015))
##d) combined
figurea7 <- grid.arrange(me_main_m2g_plot, me_main_m3d_t1_sb_plot, me_r6_m13d_plot, ncol=2, nrow=2)
ggsave(figurea7, file='../figures/figurea7.pdf', width = 14, height = 10, units="cm",dpi=1000)